Electron Self Energy for the K and L Shell at Low Nuclear Charge 
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A nonperturbative numerical evaluation of the one-photon electron self energy for the K- and 
L-shell states of hydrogenlike ions with nuclear charge numbers Z = 1 to 5 is described. Our 
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I. INTRODUCTION 

The nonperturbative numerical evaluation of radiative 
corrections to bound-state energy levels is interesting for 
two reasons. First, the recent dramatic increase in the 
accuracy of experiments that measure the transition fre- 
quencies in hydrogen and deuterium ||l|-Ql necessitates 
a numerical evaluation (nonperturbative in the binding 
Coulomb field) of the radiative corrections to the spec- 
trum of atomic systems with low nuclear charge Z . Sec- 
ond, the numerical calculation serves as an independent 
test of analytic evaluations which are based on an expan- 
sion in the binding field with an expansion parameter Za. 

In order to address both issues, a high-precision nu- 
merical evaluation of the self energy of an electron in 
the ground state in hydrogenlike ions has been per- 
formed ||,^. The approach outlined in Q is generalized 
here to the L shell, and numerical results are obtained 



and 2P3/2. Results 



Li2+, Be^+, and 



for the {n = 2) states 2Si/2, 2Pi/2 
are provided for atomic hydrogen, He^ 
B4+. 

It has been pointed out in |Q, || that the nonpertur- 
bative effects (in Za) can be large even for low nuclear 
charge and exceed the current experimental accuracy for 
atomic transitions. For example, the difference between 
the sum of the analytically evaluated terms up to the 
order of a (Zotf and the final numerical result for the 
ground state is roughly 27 kHz for atomic hydrogen and 
about 3200 kHz for He^. For the 2S state the difference 
is 3.5 kHz for atomic hydrogen and 412 kHz for He"*". 
The large difference between the result obtained by an 
expansion in Za persists even after the inclusion of a re- 
sult recently obtained in [|6| for the logarithmic term of 
order a(Za)' \Yi(Za)~'^ . For the ground state, the dif- 
ference between the all-order numerical result and the 
sum of the perturbative terms is still 13 kHz for atomic 
hydrogen and 1600 kHz for He"*". For the 2S state, the 
difference amounts to 1.6 kHz for atomic hydrogen and 
to 213 kHz for Hc+. 

These figures should be compared to the current exper- 
imental precision. The most accurately measured tran- 
sition to date is the 1S-2S frequency in hydrogen; it has 
been measured with a relative uncertainty of 1.8 parts in 
10^** or 46 Hz p|. This experimental progress is due in 
part to the use of frequency chains that bridge the range 
between optical frequencies and the microwave cesium 
time standard. The uncertainty of the measurement is 
likely to be reduced by an order of magnitude in the near 
future 1^,13] • With trapped hydrogen atoms, it should 
be feasible to observe the 1S-2S frequency with an ex- 
perimental linewidth that approaches the 1.3 Hz natural 
width of the 2S level |,§. 

The perturbation series in Za is slowly convergent. 
The all-order numerical calculation presented in this pa- 
per essentially eliminates the uncertainty from unevalu- 
ated higher-order analytic terms, and we obtain results 
for the self-energy remainder function Gse with a preci- 
sion of roughly 0.8 x Z"* Hz for the ground state of atomic 
hydrogen and 1.0 x Z'^Hz for the 2S state. 



In the evaluation, we take advantage of resummation 
and convergence acceleration techniques. The resumma- 
tion techniques provide an efhcient method of evaluation 
of the Dirac-Coulomb Green function to a relative un- 
certainty of 10^^'' over a wide parameter range 0. The 
convergence acceleration techniques remove the princi- 
pal numerical difRculties associated with the singularity 
of the relativistic propagators for nearly equal radial ar- 
guments 0. 

The one-photon self energy treated in the current in- 
vestigation is about two orders of magnitude larger than 
the other contributions to the Lamb shift in atomic hy- 
drogen. A comprehensive review of the various contri- 
butions to the Lamb shift in hydrogenlike atoms in the 
full range of nuclear charge numbers Z = 1-110 has been 
given in (ll]-|l|] . 

This paper is organized as follows. The method of eval- 
uation is discussed in Sec. O. The calculation is divided 
into a low-energy part and a high-energy contribution. 
The low-energy part is treated in Sec. HI, and the high- 
energy part is discussed in Sec. |rv|. Numerical results are 
compiled in Sec. M. Also in Sec. [V|, we compare numeri- 
cal and analytic results for the Lamb shift in the region of 
low nuclear charge numbers. Of special importance is the 
consistency check with available analytic results |[l5|,[l6| 
for higher-order binding corrections to the Lamb shift. 
We make concluding remarks in Sec. VI. 



■^-At^xinlj) In(Za) ^ + G-s,^{nlj,Zoi)\ 



(2.2) 



II. METHOD OF EVALUATION 

A. Status of Analytic Calculations 

The (real part of the) energy shift Ai^sE due to the 
electron self-energy radiative correction is usually writ- 
ten as 



AKc 



SE 



a (Za) 



■ F{nlj, Za) rric c 



(2.1) 



where i^ is a dimensionless quantity. In the following, 
the natural unit system with h — c — rua — 1 and 
e^ = Ana is employed. Note that F{nlj, Za) is a dimen- 
sionless function which depends for a given atomic state 
with quantum numbers n, I and j on only one argument 
(the coupling Za). For excited states, the (nonvanishing) 
imaginary part of the self energy is proportional to the 
(spontaneous) decay width of the state. We will denote 
here the reaZ part of the self energy by Ai5gE, exclusively. 
The semi-analytic expansion of F{nlj, Za) about Za — 
for a general atomic state with quantum numbers n, / and 
J gives rise to the semi-analytic expansion, 

F{nlj,Za) = Auinlj) ln{Za)~'^ 
+A4a{nlj) + {Za) A^oinlj) 
+ {Za)^ [Ae2{nlj) \n^{Za)-^ 



For particular states, some of the coefficients may van- 
ish. Notably, this is the case for P states, which are less 
singular than S states at the origin [see Eq. ( |2.4| ) below]. 
For the nSi/2 state {I — 0, j — 1/2), none of the terms 
in Eq. ( ^.2| ) vanishes, and we have, 

F{nSi/2,Za) = A4i(nSi/2) In(Za)"^ 
-|-A4o(nSi/2) + (Za) A5o(nSi/2) 
+ (Za)2 [A62(nSi/2) In'(^a)"' 

+A6i(nSi/2) HZa)-^ + GsE{nS,/2,Za)] . (2.3) 

The A coefficients have two indices, the first of which de- 
notes the power of Za [including those powers implicitly 
contained in Eq. (pTJ)], while the second index denotes 
the power of the logarithm ln{Za)~'^. For P states, the 
coefficients A41, A^q and Aq2 vanish, and we have 

F{nPj,Za) ^ Aio{nPj) 

+ {Za)^ [AeiinP,)\niZa)-^ + GsE(.nP,,Za)] . (2.4) 

For S states, the self-energy remainder function Gse can 
be expanded semi-analytically as 

GsE{nSi/2, Za) = Aeo{nSi/2) 

+ {Za) [A7i(nSi/2) ln(Za)-2 

+A7n{nSi/2)+o{Za)] (2.5) 

(for the "order" symbols o and O we follow the usual con- 
vention, see e.g. [E7,18|). For P states, the semi-analytic 
expansion of Gse reads 



GsE{nP,,Za)^Aeo{nPj) 
+ {Za) [yl7o(nPj)+o(Za)] 



(2.6) 



The fact that Aji{nPj) vanishes has been pointed out 
in [^. We list below the analytic coefficients and the 
Bethe logarithms relevant to the atomic states under in- 
vestigation. For the ground state, the coefficients A41 
and A40 were obtained in [M-Ga], the correction term 
A50 was found in |Pq-|28||, and the higher-order binding 
corrections Ae2 and Aqi were evaluated in [|l5|, pO]- 
The results are. 



^4l(lSi/2) = - 



^4o(lSi/2) = y-^lnA;o(lS) 



^5o(lSi/2)=47r 



139 

128 



1 



In 2 



A2(lSi/2) = -1 

28 

Al(lSl/2) = ^ 



ln2- 



21 
20 



(2.7) 



The Bethe logarithm Infco(lS) has been evaluated in ||3^ 
and igMiall as 



(2.8) 



Infco(lS) = 2.984 128 555 8(3). 
For the 2S state, we have 

^4l(2Si/2) = ^, 

A4o(2Si/2) = y-^lnA:o(2S), 



^5o(2Si/2)=47r 

A2(2Si/2) = -l 



139 

128 



1 



In 2 



Ai(2Si/2) = yln2 + ^. 



(2.9) 



The Bethe logarithm lnfco(2S) has been evaluated 
(see [p^-H, the results exhibit varying accuracy) as 



ln/co(2S) = 2.811 769 893(3). 



(2.10) 



It might be worth noting that the value for lnA;o(2S) 
given in p4| evidently contains a typographical error. 
Our in depen dent re-evaluation confirms the result given 
in Eq. ( ^.IC ) , which was originally obtained in pS] to the 
required precision. For the 2Pi/2 state we have 



^4o(2Pi/2) = -i-^lnfco(2P), 



Al(2Pl/2) 



103 
108 



(2.11) 



Note that a general analytic result for the logarithmic 
correction Aqi as a function of the bound state quan- 
tum numbers n, I and j can be inferred from Eq. (4.4a) 
of |]34|,|35|| upon subtraction of the vacuum polarization 
contribution implicitly contained in the quoted equation. 
The Bethe logarithm for the 2P states reads |3§,E3] 



lnfco(2P) = -0.030 016 708 9(3) , 



(2.12) 



Because the Bethe logarithm is an inherently nonrela- 
tivistic quantity, it is spin-independent and therefore in- 
dependent of the total angular momentum j for a given 
orbital angular momentum I. For the 2P3/2 state the 
analytic coefficients are 



^4o(2P3/2) = ^-^lnA:o(2P), 
29 

^6l(2P3/2)=go. 



(2.13) 



We now consider the limit of the function GsE{Za) as 
Za — > 0. The higher-order terms in the potential expan- 
sion (see Fig. || below) and relativistic corrections to the 
wavefunction both generat e te rms of higher order in Za 
which are manifest in Eq. ( |2.2| ) in the form of the nonva- 
nishing function GsEiZa) which summarizes the effects 
of the relativistic corrections to the bound electron wave 
function and of higher-order terms in the potential ex- 
pansion. For very soft virtual photons, the potential ex- 
pansion fails and generates an infrared divergence which 
is cut off by the atomic momentum scale, Za. This cut- 
off for the infrared divergence is one of the mechanisms 
which lead to the logarithmic terms in Eq. ( p.2| ). Some 
of t he n onlogarithmic terms of relative order {Za}'^ in 
Eq. (2^) are generated by the relativistic corrections to 
the wave function. The function Gse does not vanish, 
but approaches a constant in the limit Za -^ 0. This 
constant can be determined by analytic or semi-analytic 
calculations; it is referred to as the Aqq coefficient, i.e. 



AGo{nlj) = GsEinlj,0) . 



(2.14) 



The evaluation of the coefficient yl6o(lSi/2) has been his- 
torically problematic |la,p4-p^. For the 2S state, there 
is currently only one precise analytic result available, 

Ao(2Si/2) - -31.84047(1) [15]. (2.15) 

For the 2Pi/2 state, the analytically obtained result is 

Ao(2Pi/2) = -0.998 91(1) [16], (2.16) 



and for the 2P3/2 state, we have 



Ao(2P3/2) = -0.503 37(1) [16] 



(2.17) 



The analytic evaluations essentially rely on an expansion 
of the relativistic Dirac-Coulomb propagator in powers of 
the binding field, i.e. in powers of Coulomb interactions 
of the electron with the nucleus. In numerical evalua- 
tions, the binding field is treated nonperturbatively, and 
no expansion is performed. 



B. Formulation of the Numerical Problem 

Numerical cancellations are severe for small nuclear 
charges. In order to understand the origin of the numer- 
ical cancellations it is necessary to consider the renor- 
malization of the self energy. The renormalization pro- 
cedure postulates that the self energy is essentially the 
effect on the bound electron due to the self interaction 
with its own radiation field, minus the same effect on a 
free electron which is absorbed in the mass of the elec- 
tron and therefore not observable. The self energy of 
the bound electron is the residual effect obtained after 
the subtraction of two large quantities. Terms associ- 
ated with renormalization counterterms are of order 1 in 



the Za-expansion, w her eas the residual effect is of or- 
der (Zaf" [see Eq. (2J)]. This corresponds to a loss 
of roughly 9 significant digits a,t Z = 1. Consequently, 
even the precise evaluation of the one-photon self energy 
in a Coulomb field presented in Bffl extends only down 
to Z = 5. Among the self-energy corrections in one- 
loop and higher-loop order, numerical cancellations in 
absolute terms are most severe for the one-loop problem 
because of the large size of the effect of the one-loop self- 
energy correction on the spectrum. 
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FIG. 1. Integration contour C for the integration over the 
energy uj = E„ — 2 of the virtual photon. The contour C 
consists of the low-energy contour Cl and the high-energy 
contour Cu- Lines shown displaced directly below and above 
the real axis denote branch cuts from the photon and elec- 
tron propagator. Crosses denote poles originating from the 
discrete spectrum of the electron propagator. The contour 
used in this work corresponds to the one used in ^ . 



For our high-precision numerical evaluation, we start 
from the regularized and renormalized expression for the 
one-loop self energy of a bound electron, 

Ai?.sR= lim He' Re I ^ lS^D,,{k\K) 



A^c 



Cf 



27r / f27r)3 ^^ 



X ill) 



lim 



1 



^- /(f - 1 - "f^V^ 
' duj f dPk 



ip ) — Am 



-ie^Re / -— 



27r / f27r)3 ^' 



^M-(A:^A) 



(V' I a^' e''=-^ G{En - cu) a" e'''''" \ ^) - Am I , (2.18) 



where G denotes the Dirac-Coulomb propagator, 

1 



G{z) 



a-p + P + V ~ z 



(2.19) 



and Am is the A-dependent (cutoff-dependent) one-loop 
mass-counter term. 



Am, 



a 

■K 



InA^ 



(/3) 



(2.20) 



The photon propagator Z?^^ (fc , A) in Eq. ( 2.1S| ) in Feyn- 
man gauge reads 



D^,{e,K)^- 



gf^f 



9i^u 



fc2 - A2 + : 



(2.21) 



The contour Cp in Eq. (2.18) is the Feynman contour, 
whereas the contour C is depicted in Fig. P. The contour 
C is employed for the w-integration in the current evalu- 
ation [see the last line of Eq. ( 2.18| )]. The energy variable 
z in Eq. (2.19) therefore assumes the value 



Z = £"„ - LJ , 



(2.22) 



where En is the Dirac energy of the atomic state, and w 
denotes the complex- valued energy of the virtual photon. 
It is understood that the limit A ^ oo is taken after all 



integrals in Eq. (2.18) are evaluated 
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FIG. 2. Separation of the low-energy contour Cl into 
the infrared part Cir and the middle-energy part Cm- As 
in Fig. hi, the lines directly above and below the real axis 
denote branch cuts from the photon and electron propaga- 
tor. Strictly speaking, the figure is valid only for the ground 
state. For excited states, some of the crosses, which denote 
poles originating from the discrete spectrum of the electron 
propagator, are positioned to the right of the line Reoj = 0. 
These poles are subtracted in the numerical evaluation. 



The integration contour for the complex- valued energy 
of the virtual photon uj in this calculation is the contour 
C employed in |4^-^] and depicted in Fig. |l]. The inte- 
grations along the low-energy contour Cl and the high- 
energy contour Ch in Fig. |l| give rise to the low- and the 
high-energy contributions AEi^ and ASh to the self en- 
ergy, respectively. Here, we employ a further separation 
of the low-energy integration contour Cl into an infrared 
contour Cir and a middle-energy contour Cm shown in 
Fig. 0. This separation gives rise to a separation of the 
low-energy part ASl into the infrared part ASir and 
the middle-energy part AEm, 



AEl = AEiB. + AE] 



M ■ 



(2.23) 



For the low-Z systems discussed here, ah comphcations 
which arise for excited states due to the decay into the 
ground state are relevant only for the infrared part. Ex- 
cept for the further separation into the infrared and the 
middle-energy part, the same basic formulation of the 
self-energy problem as in ^^ is used. This leads to the 
following separation: 



cue {^En,En)±iS 

UJ £ En + i ( — 00, +00) 



infrared part AE'ir, 
middle-energy part A_Em, 
high-energy part A_Eh- 



Integration along these contours gives rise to the infrared, 
the middle-energy, and the high-energy contributions to 
the energy shift. For all of these contributions, lower- 
order terms are subtracted in order to obtain the contri- 
bution to the self energy of order (Za)'^. We obtain for 
the infrared part. 



AEi 



IR 



a 



21 



iP) 



43 



200 " ' 600 



3 Fm{nlj,Za) 



(2.24) 



where Fi^{nlj, Za) is a dimensionless function of order 
one. The middle-energy part is recovered as 



A^M = - 



200<^> + 200<^^ 

+ ^^FM{nl„Za) 



(2.25) 



and the high-energy part reads |^,Q 



A^H = Am 



|m - liv) 



+ ^—^Fii{nL,Za) 



(2.26) 



The infrared part is discussed in Sec. Ill A. The middle- 



energy part is divided into a middle-energy subtraction 
term Fma and a middle-energy remainde r Fm b- The 
subtraction term Fma is discussed in Sec. [II B| , the re- 
mainder term -Fmb is treated in Sec. IIIC. We recover 
the middle-energy term as the sum 

FM{nlj,Za) — FMA{nlj,Za) + FMB{nlj,Za) . (2.27) 

A similar separation is employed for the high-energy part. 
The high-energy part is divide d int o a subtraction term 
Fha, which is evaluated in Sec. IV A , and the hi gh-energy 
remainder Fhb , which is discussed in Sec. IV B , The sum 



of the subtraction term and the remainder is 

Fiiinlj,Za) = FnAinlj,Za)+FnB{nlj,Za). (2.28) 
The total energy shift is given as 



AE'sE — A_EiR + A_Em 
a [Za'^ 



Ey\ - Am 



[Fm{nlj,Za) 
+FMinlj,Za) + Fiiinlj, Za)] . (2.29) 



The scaled self-energy function F defined in Eq. ( |2.l| ) is 
therefore obtained as 

F(nlj, Za) — Fm(nlj, Za) 

+FMinlj,Za) + Fiiinlj,Za). (2.30) 

In analogy to the approach described in |^, ^^, we 
define the low-energy part as the sum of the infrared part 
and the middle-energy part. 



ASl = AEiR + AEi 



M 



l{f3) + l{V)+^-^FUnl„Za) 



, (2.31) 



where 



Fi^inlj, Za) = Fm{nlj,Za) + FM{nlj,Za) . (2.32) 



The limits for the functions Fi^(nL, Za) and Fn^nlj, Za) 
as Za — > were obtained in [5 , 48 , 50| . 



C. Treatment of the divergent terms 

The free electron propagator 
F 



a ■ p + (3 — z 



(2.33) 



and the full electron propagator G defined in Eq. (2.19) 
fulfill the following identity, which is of particular impor- 
tance for the validity of the method used in the numerical 
evaluation of the all-order binding correction to the Lamb 
shift: 



G = F- FVF + FVGVF. 



(2.34) 



This identity leads naturally to a separation of the one- 
photon self energy into a zero- vertex, a single- vertex, and 
a many-vertex term. This is represented diagrammati- 
cally in Fig. || 




^ ^ 



FIG. 3. The exact expansion of the bound electron prop- 
agator in powers of the binding field leads to a zero-potential, 
a one-potential, and a many-potential term. The dashed lines 
denote Coulomb photons, the crosses denote the interaction 
with the (external) binding field. 

All ultraviolet divergences which occur in the one- 
photon problem (mass counter term and vertex diver- 
gence) are generated by the zero-vertex and the single- 
vertex terms. The many-vertex term is ultraviolet safe. 
Of crucial importance is the observation that one may 
additionally simplify the problem by replacing the one- 
potential term with an approximate expression in which 
the potential is "commuted to the outside." The ap- 
proximate expression generates all divergences and all 
terms of lower order than a (Za)^ present in the one- 
vertex term. Unlike the raw one-potential term, it is 
amenable to significant further simplification and can be 
reduced to one-dimensional numerical integrals that can 
be evaluated easily (a straightforward formulation of the 
self-energy problem requires a i/iree-dimensional numer- 
ical integration). Without this significant improvement, 
an all-order calculation would be much more difficult at 
low nuclear charge, because the lower-order terms would 
introduce significant further numerical cancellations. 

In addition, the special approximate resolvent can be 
used effectively for an efficient subtraction scheme in the 
middle-energy part of the calculation. In the infrared 
part, such a subtraction is not used because it would in- 
troduce infrared divergences. 

We now turn to the construction of the special ap- 
proximate resolvent, which will be referred to as Ga and 
will be used in this calculation to isolate the ultraviolet 
divergences in the high-energy part (and to provide sub- 
traction terms in the middle-energy part). It is based on 
an approxi matio n to the first two terms on the right-hand 
side of Eg. (|2.34 ) . The so-called one-potential term FVF 
in Eq. (2.34) is approximated by an expression in which 
the potential terms V are commuted to the outside: 



~FVF^~^ {V,F^} . 
Furthermore, the following identity is used: 



(2.35) 



F^ = 



1 



a. ■ p + P — z ^ 

1 2z{l3 + z) 



p^ + l-Z^ (p2 + 1 _ ^2 



+ - 



2z{a ■ p) 
(p2 + 1 _ z^' 



(2.36) 



In 2 X 2 spinor space, this expression may be divided into 
a diagonal and a non-diagonal part. The diagonal part 
is 



diag(F2) = 



1 



2 z (/? + z) 



p2 + 1 - z2 {p^ + 1- z2) 

The off-diagonal part is given by 

2 z {a ■ p) 



(2.37) 



F^ - diag(i^^) 



(p2 + 1 _ ^2)2 



We define the resolvent Ga as 



1 



Ga=F-- {F,diag(^2)| 



(2.38) 



All divergences which occur in the self energy are gen- 
erated by the simplified propagator Ga. We define the 
propagator Gb as the difference of G and Ga , 

Gb = G — Ga 

= - {V,diag{F^)} -FVF + FVGVF. (2.39) 

Gb does not generate any divergences and leads to the 
middle-energy remainder discussed in Sec. Ill C| and the 
high-energy remainder (Sec. IV B|). 



III. THE LOW-ENERGY PART 
A. The Infrared Part 

The infrared part is given by 

X (i^ I a'^ e"=-" G{E^ - u) a" e^""'" | ^) , 



(3.1) 



where relevant definitions of the symbols can be found in 
Eqs. (2.18-2.21), the contour Gir is as shown in Fig. H, 
and the unregularizcd version of the photon propagator 



D^^Ak^) 



9t^i^ 



e 



(3.2) 



may be used. The infrared part consists of the following 
integration region for the virtual photon: 



uj€ (0,^ E„) ±iS ] 

; „ ^° \ } infrared part AE^ ■ (3.3) 

Following Sees. 2 and 3 of [O, we write AE'ir as a 
three-dimensional integral [see, e.g., Eqs. (3.4), (3.11), 
and (3.14) of @] 

Ai^iH = ^ % - ^ (P.V.) /"" d. 



dxi a;;^ / dx2 a;2 A4ir(x2,xi, z) , (3.4) 

^0 



is made, so that r G (0,1) and y € (0, oo). The scaling 
variable a is defined as 



The Jacobian is 

i9(a;2,xi) 



d{r,y) 



271^^2. 



dx2 dxi 

dr dr 

dx2 dxi 

~dy ~dy 



The function S'ir is given by 
2r'^y^ 



where 



SiK{r,y,z) = -- 



Mm [ — , -, 
\ a a 



(3.9) 



(3.10) 



Mm,{x2,xi,z) =X! X! 

MX2) G'^{X2,X1,Z) fj{xi) A'^{X2,XI) . (3.5) 

Here, the quantum number k is the Dirac angular quan- 
tum number of the intermediate state, 



^ = 2{l-j) (j + 1/2) , 



(3.6) 



where I is the orbital angular momentum quantum num- 
ber and j is the total angular momentum of the bound 
electron. The functions fi{x2) {i = 1 , 2) are the radial 
wave functions defined in Eq. (A. 4) in |47|| for an arbi- 
trary bound state (and in Eq. (A. 8) in ||47t for the IS 
state). We define i = 3 — i. The functions G]^{x2,xi, z) 
{i,j = 1,2) are the radial Green functions, which re- 
sult from a dec ompo sition of the electron Green function 
defined in Eq. ( ^.19 ) into partial waves. The explicit for- 
mulas are given in Eq. (A. 16) in |l7| . 

The photon angular functions A]^ {i,j = 1,2) are de- 
fined in Eq. (3.15) of Rcf. |^^ for an arbitrary bound 
state. In Eq. (3.17) in Q, specific formulas are given 
for the IS state. In Eqs. (2.2), (2.3) and (2.4) of [g, 
the special cases of S1/2, P1/2 and P3/2 states are con- 
sidered. Further relevant formulas for excited states can 
be found in |5l| ]. The photon angular functions depend 
on the energy argument z, but this depende nce i s usually 
suppressed. The summation over k in Eq. (3^) extends 
over all negative and all positive integers, excluding zero. 
We observe that the integral is symmetric under the in- 
terchange of the radial coordinates X2 and xi , so that 



. „ a En 2 a , 

AEm = - ^ P.V. 

TT 10 n 



dz 



10 



E„ 



/ dxixl dx2xl Mm{x2,xi,z) . (3.7) 
Jo Jo 

The following variable substitution. 



(.. o r 00 2 

^-^E E E/.(?) 

\ a a / \a/ \ a a/ 

2^2 5 °o 
= ^— '^Tij^,\K\ir,y,z), (3.11) 

|k.| = 1 

where in the last line we define implicitly the terms Tjj^ 1^1 
for |k| = 1, . . . , 00 as 



Tm,\^,\{r,y,z)= ^ ^ 

k=±|k;| i,j = l 

h(^)G^^Az)fJy)A^(:i,y). (3.12) 



Using the definition (3.11), we obtain for AEiji 



a E„ a , , /" " 
Ai^iR - - -^ + - P.V. / dz 

TT 10 TT J-^E^ 



dr I dySm(r,y,z) 
Jo 



(3.13) 



The specification of the principal value (P.V.) is neces- 
sary for the excited states of the L shell, because of the 
poles along the integration contour which correspond to 
the spontaneous decay into the ground state. Here we 
are exclusively concerned with the real part of the energy 
shift, as specified in Eq. (|3.l|) , which is equ ivalent to the 
specification of the principal value in (3.13). Evaluation 
of the integral over z is facilitated by the subtraction of 
those terms which generate the singularities along the in- 
tegration contour (for higher excited states, there can be 



numerous bound state poles, as pointed out in |51 



X2IX1, y^axi, 



(3.8) 



For the 2S and 2Pi/2 states, only the pole contribution 
from the ground state must be subtracted. For the 2P3/2 
state, pole contributions originating from the IS, the 2S 
and the 2Pi/2 states must be taken into account. The nu- 
merical evaluation of the subtracted integrand proceeds 



along ideas outlined in ^,Q and is not discussed here 
in any further detail. 

The scaling parameter a for the integration over y is 
chosen to simplify the exp onen tial dependence of the 
function S defined in Eq. (3.11). The main exponen- 
tial dependence is given by the relativistic radial wave 
functions (upper and lower components). Both compo- 
nents [/i(a;) and f2ix)] vary approximately as (neglecting 
relatively slowly varying factors) 



exp {—ax/2) 



(for large x) . 



The scaling variable a, expanded in powers of Za, is 



a = 2^l-E^ 



{Zaf 
2v? 



O [{ZaY 



= 2— + 0[{Zaf] 
n 



(3.14) 



Therefore, a is just twice the inverse of the Bohr radius 
n/{Za) in the nonrelativistic limit. The product 



!<i)^ 



f'il) 



for arbitrary i, j € {1, 2} 



[which occurs in Eq. (3.11)] depends on the radial argu- 
ments approximately as 



X exp 



(l-r)y 



(for large y) . 



Note that the main dependence as given by the term 
exp(— y) is exactly the weight factor of the Gaufi- 
Laguerre integration quadrature formula. The deviation 
from the exact exp(— j/)-type behavior becomes smaller 
as r ^ 1. This is favorable because the region near r = 1 
gives a large contribution to the integral in (3.13). 



TABLE I. Infrared part for the K and L shell states, Fir(1Si/2, ^a), Fin{2Si/2, Za), Fin{2Pi/2, Za), and 
FiR{2P3/2,Za), evaluated for \ow-Z hydrogenlike ions. The calculations were performed with the numerical value 
of Q~^ = 137.036 for the fine-structure constant. 



z 


FiR(lSi/2,^a) 


Fm{2Si,2,Za) 


Fm{2Pi/2,Za) 


Fm(2P3/2,^Q) 


1 

2 
3 

4 
5 


7.236 623 736 8(1) 
5.539 002 119 1(1) 
4.598 155 821 8(1) 
3.963 124 140 6(1) 
3.493 253 319 4(1) 


7.479 764 180(1) 
5.782 025 637(1) 
4.840 923 962(1) 
4.205 501 798(1) 
3.735 114 958(1) 


0.085 327 852(1) 
0.086 073 669(1) 
0.087 162 510(1) 
0.088 543 188(1) 
0.090 180 835(1) 


0.082 736 497(1) 
0.083 279 461(1) 
0.084 091 830(1) 
0.085 140 788(1) 
0.086 403 178(1) 



The sum over |k| in Eq. ( |3.11 ) is carried out locally, 
i.e., for each set of arguments r, y, z. The sum over \k\ is 
absolutely convergent. For \k\ — > cxd, the convergence of 
the sum is governed by the asymptotic behavior of the 
Bessel functions which occur in the photon functions A^ 
{i,j = 1,2) [see Eqs. (3.15) and (3.16) in [Q]- The pho- 
ton functions contain products of two Bessel functions of 
the form J7i(p2/i) where J7i stands for either ji or j^', and 
the index I is in the range I € {|k| — 1, |k|, |k| + 1}. The 
argument is either p2 = [En — z) X2 or pi = [En — z) xi. 
The asymptotic behavior of the two relevant Bessel func- 
tions for large I (and therefore large |k|) is 



fii^) = - 



I a;' 



Ji{x) 



X (2^ + 1)!! 
(21 + 1)11 



1 + 0|1 



and 



1 + 



(3.15) 
(3.16) 



This implies that when min{p2,Pi} — P2 < I, the func- 
tion Si{p2) vanishes with increasing I approximately as 
{ep2/2Vf. This rapidly converging asymptotic behavior 
sets i n as soon as I k, \k\ > p2 = rtoy/a [see Eqs. (2.22) 
and ( |3.12 )]. Due to the rapid convergence for \k\ > p2, 



the maximum angular momentum quantum number \k,\ 
in the numerical calculation of the infrared part is less 
than 3 000. Note that because z E {-^ En,En) in the 
infrared part, uj < jq En. 

The integration scheme is based on a crude estimate 
of the depe ndence of the integrand SiB,{r,y, z) defined 
in Eq. (3.11) on the integration variables r, y and z. The 
main contribution to the integral is given by the region 
where the arguments of the Whittaker functions as they 
occur in the Green function [see Eq. (A. 16) in |4^] are 
much larger than the Dirac angular momentum, 

2c^ > \k\ 
a 

(see also p. 56 of [M). We assume the asymptotic form 
of the Green function given in Eq. (A. 3) in [Q applies 
and attribute a factor 



exp[-(l ~r)cy/a] 



radial Green functions GJf as they occur in 

Note that relatively slowly varying factors 

The products of the radial wave 



to the 

Eq. ( pl| ) 

are replaced by unity. 



func tions /j and /j, according to the discussion following 
Eq. ( 3.14 ), behave as 



exp 



{l-r)y 



for large y. The photon functions A]^ in Eq. ( 3.11 ) give 
rise to an approximate factor 



sin[(l - r) {En ~ z) y/a] 
(1-r) 



(3.17) 



Therefore [see also Ec^. (2.12) in [Q], we base our choice 
of the integration routine on the approximation 



exp 



c 1 



(1 -r)j/ 



sin[(l -r) [En - z)y/a\ 
(1-r) 



(3.18) 



for S'iR. The three-dimensional integral in ( 3.13 ) is eval- 
uated by successive Gaussian quadrature. Details of the 
integration procedure can be found in ||5|]. 

In order to check the numerical stability of the results, 
the calculations are repeated with three different values 
of the fine-structure constant a: 



a< = 1/137.036 000 5, 
ao = 1/137.036 000 
a> 



1/137.035 999 5. 



and, 



(3.19) 



These values are close to the 1998 CODATA recom- 
mended value of a-i = 137.035 999 76(50) ||. The 
calculation was parallelized using the Message Passing 
Interface (MPI) and carried out on a cluster of Silicon 
Graphics workstations and on an IBM 9276 SP/2 multi- 
processor system [ p4[ . T he results for the infrared part 
Fir, defined in Eq. ( 2.24| ), are given in Table for a value 
of a~^ = a^^ = 137.036. This value of a will be used 
exclusively in the numerical evaluations presented here. 
For numerical results obtained by employing the values 
of a< and a> [see Eq. ( 3.19| )] we refer to [||. 



B. The Middle-Energy Subtraction Term 



The middle-energy part is given by 



AEm^ 



Cm 



2n / (27r) 



D^.ik') 



X (V' I a" e' '=■"' GiEn - lo) a'' e~' '=-=^ | V) , (3.20) 

where relev a nt de finitions of the symbols can be found in 
Eqs. (|2.18|)-(p^ and Eq. (jsj), and the contour Cm is 
as shown in Fig. 0. The middle-energy part consists of 
the following integration region for the virtual photon: 



^G (Tr;En,En)±iS 



e {0,^En)±iS 



middle-energy part AE] 



M ■ 



(3.21) 



The numerical evaluation of the middle-energy part is 
simplified considerably by the decomposition of the rcla- 
tivistic Dirac-Coulomb Green function G as 



G 



Ga + Gb , 



(3.22) 



where Ga is defined in (2.38) and represents the sum of 
an approximation to the so-called zero- and one-potential 
terms generated by the expansion of the Dirac-Coulomb 
Green function G in powers of the binding field V. We 
define the middle-energy subtraction term Fma as the 
expression obtained upon substitution of the propagator 
Ga for G in Eq. ( 3.20| ). The propagator Gb is simply 
calculated as the difference of G and Ga [see Eq. ( 2.39 )]. 
A substitution of the propagator Gb for G in Eq. ( 3.20] ) 
leads to the middle -energy remainder Fmb which is dis- 
cussed in Sec. IIIC| . We provide here the explicit expres- 



9 / doj f d'^k „ ,,9, 



^3 ^p.u 



(V- 



f-i i k-x 



a^'e 



Cm 27r J (2^) 



GA(F„-w)a'^e-"=-=^|V) (3.23) 



and 



AE, 



MB — — le 



dw f d^fe 



D^,,{e 



Cm 2^ J (2^)' 

X (1/; I a" e"=-=" Gb(F„ - uj) a" e"""-^ | V) ■ (3.24) 

Note that the decomposition of the Dirac-Coulomb Green 
function as in ( 3.22D is not applicable in the infrared part, 
because of numerical problems for ultra-soft photons (in- 
frared divergences) . Rewriting (|3.23D appropriately into 
a three-dimensional integral [H, [47|, [48|| , we have 



AE- 



a 9 



MA 



TT 10 TT 



2 a /"lo 



dz 



/■oo nxi 

/ dxix^ / dx2xl Mma{x2,Xi,z) . 
Jo Jo 



(3.25) 



The function A^ma(2;2j 2:1, z) is defined i n an alogy to the 
function A4i^{x2,xi, z) defined in Eq. ( |3.5| ) for the in- 
frared part. Also, we define a function SMAix2,xi, z) in 



analogy to the function 5iR(a:2: Xi, z) given in Eq. ( 3.11 ) 
for the infrared part, which will be used in Eq. (3.25) 
below. We have. 



SuAir, y, z) = 



2 r2 2/5 



Mm A 



\ a a / 



|k.| = 1 k;=±|k;| ij = l 



With these definitions, the middle-energy subtraction 
term Ai?MA can be written as 



^£,1 



A^MA - - 77: £^n + - / dz 

TT ID TT 



dy / dr S'MA(r, zj, z) . (3.28) 

Jo 



(3.26) 



2 ^2 5 "^ 
a° ^ — ' 

|k| = 1 

The expansion of the propagator Ga into partial waves 



The subtracted lower-order terms yield, 



is 



and in Eqs. (D.37) 



given in Eqs. (5.4) and (A. 20) in [^ and in Eqs. (D.37 
and (D.42) in H. This expansion leads to the compo- 
t functions GY „. The terms Tma i^i in the last line 



A^MA = - 
TT 



279,^, 219 ,_,^, 

200<^> + 200<^> 



nent functions G)^^. The terms Tma.Ik] in 
of Eq. ( |3^ ) read ' 



-F 



(Za)4 



FMA{nlj,Za) 



(3.29) 



k.=±|k| ij' = l 



A 



(?)<=-(?■!-) MD-^e 



i^. M) . (3.27) 
" aJ 



The three-dimensional integral in (3.28) is evaluated by 
successive Gaussian quadrature. Details of the integra- 
tion procedure can be found in p| . The numerical results 
are summarized in Table O. 



TABLE II. Numerical results for the middle-energy subtraction term -Fma, the middle-energy 
remainder term -Fmb, and the middle-energy term Fm- The mid dle-en er gy te rm Fm is gi ven as the 
sum Fm{nlj,Za) — FMA{nlj,Za) + FMB{nlj, Za) [see also Eqs.(p.25|), (p.29|), and (p.33|)]. 



z 


-Fma(1Si/2,^q) 


-fMA(2Si/2, Za) 


Fma{2Pi/2, Za) 


FMA{2P3/2,Za) 


1 

2 
3 

4 

5 


2.699 379 904 5(1) 
2.659 561 381 1(1) 
2.623 779 453 0(1) 
2.591 151 010 1(1) 
2.561 096 522 1(1) 


2.720 878 318(1) 
2.681 820 660(1) 
2.647 262 568(1) 
2.616 290 432(1) 
2.588 297 638(1) 


0.083 207 314(1) 
0.084 208 832(1) 
0.085 831 658(1) 
0.088 040 763(1) 
0.090 803 408(1) 


0.701 705 240(1) 
0.701 850 024(1) 
0.702 091 147(1) 
0.702 426 850(1) 
0.702 854 461(1) 


Z 


-Fmb(1Si/2,^q) 


Fmb(2Si/2,^q) 


FMBi2Pi/2, Za) 


Fmb{2P3/2, Za) 


1 
2 
3 

4 

5 


1.685 993 923 2(1) 
1.626 842 294 5(1) 
1.571 406 090 7(1) 
1.519 082 768 6(1) 
1.469 482 409 0(1) 


1.784 756 705(2) 
1.725 583 798(2) 
1.670 086 996(2) 
1.617 650 004(2) 
1.567 873 140(2) 


0.771 787 771(2) 
0.770 778 394(2) 
0.769 153 314(2) 
0.766 954 435(2) 
0.764 220 149(2) 


-0.094 272 681(2) 
-0.094 612 071(2) 
-0.095 165 248(2) 
-0.095 922 506(2) 
-0.096 874 556(2) 


Z 


FM(lSi/2,^a) 


FM(2Si/2,^a) 


FM{2Pi/2,Za) 


FM{2P3/2,Za) 


1 
2 
3 
4 
5 


4.385 373 827 7(1) 
4.286 403 675 7(1) 
4.195 185 543 6(1) 
4.110 233 778 8(1) 
4.030 578 931 1(1) 


4.505 635 023(2) 
4.407 404 458(2) 
4.317 349 564(2) 
4.233 940 436(2) 
4.156 170 778(2) 


0.854 995 085(2) 
0.854 987 226(2) 
0.854 984 972(2) 
0.854 995 198(2) 
0.855 023 557(2) 


0.607 432 559(2) 
0.607 237 953(2) 
0.606 925 899(2) 
0.606 504 344(2) 
0.605 979 905(2) 



C. The Middle-Energy Remainder 

The remainder term in the middle-energy part involves 
the propagator Gb defi ned i n Eq. ( |2.39D , Gb = G — Ga , 
where G is defined in (2.19) and Ga is given in (2.38). 



In analogy to the middle-energy subtraction term, the 
middle-energy remainder can be rewritten as a three- 
dimensional integral. 



10 



AS, 



9 p 



MB 



/•l /"OO 

dr dySMB{r,y,z), (3.30) 



where 



2 ,,5 



SMB{r,y,z) = — ^ E E E /^(v) 

The functions Gg ^ are obtained as the difference of the 
expansion of the full propagator G and the simplified 
propagator Ga into angular momenta, 






(3.32) 



where the G^ are listed in Eq. (A. 16) in |^^ and in 
Eq. (D.43) in |^, and the G^ ^ have already been de- 
fined in Eqs. (5.4) and (A.20) in ^ and in Eqs. (D.37) 
and (D.42) in [g|. There are no lower-order terms to sub- 
tract, and therefore 



AE- 



(Za) 



MB 



■ FMB{nlj,Za) 



(3.33) 



The three-dimensional integral ( 3.30| ) is evaluated by suc- 
cessive Gaussian quadrature. Details of the integration 
procedure are provided in [^. Numerical results for the 
middle-energy remainder -Fmb are summarized in Ta- 
ble H for the K- and L-shell states. 

For the middle-energy part, the separation into a sub- 
traction and a remainder term has considerable compu- 
tational adva ntage s whi ch be come obvious upon inspec- 
tion of Eqs. ( 3.29 ) and ( 3.33 ). The subtraction involves 
a propagator whose angular components can be evalu- 
ated by recursion [plEql, which is computationally time- 
consuming. Because the subtraction term involves lower- 



order components [see Eq. (2.25)], it has to be evaluated 
to high precision numerically (in a typical case, a rela- 
tive uncertainty of 10~^^ is required). This high preci- 
sion requires in turn a large number of integration points 
for the Gaussian quadratures, which is possible only if 
the numerical evaluation of the integrand is not compu- 
tationally time-consuming. For the remainder term, no 
lower-order terms have to be subtracted, and the rela- 
tive precision required of the integrals is in the range 
of 10^^^. ..10^^. A numerical evaluation to this lower 
level of precision is feasible, although the calculation of 
the Green function Gb is computationally more time con- 
suming than that of Ga [^,47,4^. The separation of the 
high-energy part into a subtraction term and a remain- 
der term, which is discussed in Sec. [V, is motivated by 



analogous considerations as for the middle-energy part. 
In the high-energy part, this separation is even more im- 
portant than in the middle-energy part, because of the 
occurrence of infinite terms which need to be subtracted 
analytica lly before a numerical evaluation can proceed 
[seeEq. (Q below]. 

We now summarize the results for the middle-energy 
part. The middle-energy part is the sum of the middle- 
energy subtraction term F ma an d the middle-energy re- 
mainder Fmb [see also Eq. ( |2.27 )]. Numerical results are 
summarized in Table O for the K- and L-shell states. The 
low-energy part Fj^ is defined as the sum of the infrared 
contribut ion f iR and the middle-energy contribution Fm 
[see Eq. ( 2.32 )]. The results for Fl are provided in the 
Table III for the K- and L-shcU states. The limits for the 



low-energy part as a function of the bound state quantum 
numbers can be found in Eq. (7.80) of [p|: 



Fi^{nlj, Za) = - 5ifl In(Za) ^ 
-llnM«,0+(ln2-ii)i 



„, „ 16\ 1 
2 In 2 



3 r 

l^ J 21 + I ' \2 ^ ^ 1) K(2/-hl) 



5\n_2|-l 

ey n{2i + i) ^ ' 



(3.34) 



The limits for the states under investigation in this paper 
are 

Fl(1Si/2, Za) = (4/3) ln(Za)"2 _ ^ 554 g42 + 0(Za) , 
i^L(2Si/2, Za) = (4/3) ln(Za)-2 _ 1.191 497 + 0{Za) , 



Fl(2Pi/2, Za) = 0.940 023 + 0{Za) , 
Fl(2P3/2, Za) = 0.690 023 + 0{Za) . 



(3.35) 



These limits are consistent with the numerical data in 
Table III. For S states, the low-energy contribution Fl 
diverges logarithmically as Za -^ 0, whereas for P states, 
Fl approaches a constant as Za -^ 0. The leading loga- 
rithm is a consequence of an infrared divergence cut off 
by the atomic momentum scale. It is a nonrelativistic 
effect which is generated by the nonvanishing probabil- 
ity density of S waves at the origin in the nonrelativistic 
limit. The presence of the logarithmic behav ior for S 
state s [nonvanishing A4i-coeffLcient, see Eqs. ( ^■2\ ) and 
( p.3[ )] and its absence for P states is reproduced consis- 
tently by the data in Table III. 
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TABLE III. Low-energy part Fl for the K- and L-shell states Fl(1Si/2, ^a), Fl(2Si/2, ^a), 
_Fl(2Pi/2, -^o), and _Fl(2P3/2, ^a), evaluated for low-Z hydrogenlike ions. 



Fl(1Si/2,Zq) 



FL(2Si/2,^a) 



i^L(2Pi/2,^a) 



FL(2P3/2,^a) 



1 


11.621 997 564 5(1) 


11.985 399 203(2) 


0.940 322 937(2) 


0.690 169 056(2) 


2 


9.825 405 794 7(1) 


10.189 430 095(2) 


0.941 060 895(2) 


0.690 517 414(2) 


3 


8.793 341 365 4(1) 


9.158 273 526(2) 


0.942 147 482(2) 


0.691 017 729(2) 


4 


8.073 357 919 4(1) 


8.439 442 234(2) 


0.943 538 386(2) 


0.691 645 132(2) 


5 


7.523 832 250 6(1) 


7.891 285 736(2) 


0.945 204 392(2) 


0.692 383 083(2) 



IV. THE HIGH-ENERGY PART 



A. The High-Energy Subtraction Term 



The high-energy part is given by 



AEv 



lim i e / —— 



A- 



duj f (Pk 



Ch 



27r 



(V-l 



a^'e' 



'■G{E„ 



(27r)3 "''' 



-DM-(fc',A) 



where relevant definitions of the symbols can be found 
in Eqs. ( ^.18 )-(2.21), and the contour Ch is as shown 
in Fig. |l|. The high-energy part consists of the following 
integration region for the virtual photon, 



ll> € {En — ioo, En + ioo) 
z e (—ioo, ioo) 



■ high-energy part AE^ . 

(4.2) 



The separation of the high-energy part into a subtraction 
term and a remainder is ac complished as in the middle- 
energy part [see Eq. ( 3.22 )] by writing the full Dirac- 
Coulomb Green function G [Eq. (|]l|)] as G = Ga + Gb. 
We define the high-energy subtraction term Fha as the 
expression obtained upon substitution of the propagator 
Ga for G in Eq. (4.1), a nd a substitution of the prop- 
agator Gb for G in Eq. (4.1) leads to t he hi gh-energy 
remainder i^HB which is discussed in Sec. [VB| . The sub- 
traction term (including all divergent contributions) is 
generated by Ga, the high-energy remainder term corre- 
sponds to Gb. We have 



Ai?HA = — lini 



Ch 



duj 

2^ 



fk_ 



^^.(A:^A) 



X (V' I a^ e"^-^ GA{En - uj) a"" q-'^'' \ ^) 



(4.3) 



and 



X (V I a'' e"=-=^ GB{En - cj) a" c-'^'^ \ ^) 



The contribution Ai?HA corresponding to Ga can be sep- 
arated further into a term AE\^'p^, which contains all di- 

vergent contributions, and a term Ai?^^, which contains 
contributions of lower order than (Za)^, but is conver- 
gent as A ^ 00. This separation is described in detail 
in EtI |o| . We have 



r(i) 



^(2) 



A^HA = AE'^l + AE^j^ . 



?W 



(4.5) 



(4 1) ^^ obtain for AiJjj^, which contains a logarithmic di- 



vergence as A ^ 00, 



AE, 



(1) 

HA 



a 



|lnA^J)w + (iln2-iI)w 



+ 



(Za)'^ 



Fi'linl,,Za) 



(4.6) 



For the contribution F^l , an explicit analytic result is 
given in Eq. (4.15) in [E^. This contribution is therefore 

not discussed in any further detail here. The contribution 

(2) 
Ai^jj^ contains lower-order terms: 



AE, 



(2) 
HA 



a 

IT 



(Za)4 



P'ilinlj.Za) 



(4.7) 



Altogether we have 



A^HA = AE, 



(1) 

HA 



AE^ 



(2) 
HA 



a 

TT 



In A' 



(/3) - I {V) 



6 



:r—FB_A[nlj,Za) 



The scaled function i^HA(^^j, ^c^) is given by 



(4.8) 



FHA{nlj,Za) = F^^linlj, Za) + FJ^l{nlj, Za) . (4.9) 
The term Ai?^^ falls naturally into a sum of four con- 



(4.4) tributions [0 
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AK 



(2) 
HA 



T1+T2+T3 + T4, 



where 



T. = --^iV) 



{Zaf 



To 



T.= 



l-^^^^M^> 



hi{nlj, Za) 
(Zar 



{Zaf 



(4.10) The hmits of the contributions Fj^^{nlj,Za) and 
Fj^/^{nlj,Za) as (Za) -^ have been investigated 
in ||^,|9[|§. For the contribution F^l{nlj,0), the re- 
sult can be found in Eq. (3.5) in E^. The hmits of the 
functions hi{nlj, Za) (i = 1, 2, 3, 4)as Za — *■ are given 
as a function of the atomic state quantum numbers in 
h2[nlj,Za) , gq^ ^33^ ^^ [gg| p^j. ^j^g scaled high-energy subtraction 

term Fha, the limits read (see Eq. (3.9) in |{49| ) 



h^{nlj, Za) , 



Ti — — hii{nlj, Za) 



(4.11) 



The functions hi {i — 1, 2, 3, 4) are defined in Eqs. (4.18), 
(4.19) and (4.21) in |47) (see also Eq. (3.6) in Q). The 
evaluation of the high-energy subtraction term proceeds 
as outlined in p7J-p9|, albeit with an increased accu- 
racy and improved calculational methods in intermediate 
steps of the calculation in order to overcome the severe 
numerical cancellations in the \ow-Z region. We recover 
Fj^J^ as the sum 



,(2) 



F^J^{nlj, Za) = hi{nlj, Za) + h2{nlj, Za) 
+h3{nlj, Za) + hi{nlj, Za) . 



(4.12) 



The scale d function FuAinlj, Za) [see also Eqs. ( |2.26| ) 
and (2.28)] is given by 



.(1) 



(2), 



FiiAinlj, Za) = F^l{nlj,Za) + F^l{nlj,Za) . (4.13) 



FYiA{nl;j,Za) 



11 
10 



2 In 2 



-In 2 
1 



In 2 



21 
1 



1 



1 ■ V2 
+ 0{Za). 



In 2 



k{21 + 1) 
(4.14) 



Therefore, the explicit forms of the limits for the states 
under investigation in this paper are 



-FHA(lSi/2,^a) = 


-1.219 627 + 0(Za) 


FHA(2Si/2,Za) = 


-1.423 054 + 0(Za) 


FHA(2Pi/2,^a) = 


-1.081 204 + O(Za) 



Fha(2P3/2, Za) = -0.524 351 + 0{Za) . 



(4.15) 



Numerical results for Fua, which arc presented in Ta - 
ble rV, exhibit consistency with the limits in Eq. ( 4.15| ). 



TABLE IV. Numerical results for the high-energy subtraction term Fha and the high-energy remainder term -Fhb- 
The high-energy term Fh is the sum Fninlj, Za) = FnAijilj, Za) + FHB{nlj, Za). 



Fha(1Si/2, Za) 



Fha{2Si/2, Za) 



-FHA(2Pi/2,Za) 



FHAi2P 3 /2, Za) 



1 

2 
3 

4 
5 


-1.216 846 660 6(1) 
-1.214 322 536 9(1) 
-1.212 026 714 1(1) 
-1.209 942 847 4(1) 
-1.208 059 033 6(1) 


-1.420 293 291(1) 
-1.417 829 864(1) 
-1.415 635 310(1) 
-1.413 693 422(1) 
-1.411 992 480(1) 


-1.081 265 954(1) 
-1.081 451 269(1) 
-1.081 760 224(1) 
-1.082 192 995(1) 
-1.082 749 845(1) 


-0.524 359 802(1) 
-0.524 385 053(1) 
-0.524 427 051(1) 
-0.524 485 727(1) 
-0.524 561 017(1) 


Z 


FHB(lSi/2,2'a) 


FHB{2Si/2,Za) 


FuB(2P^/2,Za) 


Fhb{2P3/2, Za) 


1 
2 
3 

4 
5 


-0.088 357 254(1) 
-0.082 758 206(1) 
-0.076 811 229(1) 
-0.070 590 991(1) 
-0.064 146 139(1) 


-0.018 280 727(5)" 
-0.012 729 99(1) 
-0.006 861 02(1) 
-0.000 746 40(1) 
0.005 567 16(1) 


0.014 546 64(1) 
0.014 574 21(1) 
0.014 620 51(1) 
0.014 685 82(1) 
0.014 770 52(1) 


-0.042 310 69(1) 
-0.042 296 81(1) 
-0.042 273 58(1) 
-0.042 240 92(1) 
-0.042 198 76(1) 


Z 


FH(lSi/2,^a) 


Fn{2S,/2,Za) 


Fh(2Pi/2,2'q) 


Fii{2P3/2,Za) 


1 
2 
3 

4 
5 


-1.305 203 915(1) 
-1.297 080 743(1) 
-1.288 837 943(1) 
-1.280 533 839(1) 
-1.272 205 173(1) 


-1.438 574 018(5) 
-1.430 559 85(1) 
-1.422 496 33(1) 
-1.414 439 82(1) 
-1.406 425 32(1) 


-1.066 719 31(1) 
-1.066 877 06(1) 
-1.067 139 72(1) 
-1.067 507 18(1) 
-1.067 979 33(1) 


-0.566 670 50(1) 
-0.566 681 86(1) 
-0.566 700 63(1) 
-0.566 726 65(1) 
-0.566 759 78(1) 



'^Result obtained with a greater number of integration nodes than are used for the higher-Z results. 
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B. The High-Energy Remainder 

The remainder term in the hig h-ener gy part involves 
the propagator Gb defi ned i n Eq. ( 2.39| ), Gb = G — Ga , 
where G is defined in ( 2.19 ) and Ga is given in ( 2.38 ). 
The energy shift is 



■ 1*1 OG /"OO 

AShb = / dz / dxi x\ 

T^ Jo Jq 

/ dx2a;2 {A^HB(a;2,a;i,z) + c.c.} , (4.16) 
Jo 

where c.c. denotes the complex conjugate. The photon 
energy integration is evaluated with the aid of the sub- 
stitution 



1 n 

z — > 1 u where u = — t 

2 V t 



(4.17) 



In analogy with the middle-energy subtr action and 



remainder terms discussed in S ees. IIIB and [II C 
[see especially Eqs. (3.26) and ( p. 31 )], the functions 



■A4hb(2:2, xi, z) and S'hb(?', y, z) and the terms Thb.IkI are 
defined implicitly in the following: 



Smi{r,y.,t) 



h) 


^2 5 

— |-Re 


1^ 


r y' V 



M 



HB 



(11. 

\ a '■ 



i2 



E E E 

|k| = 1 k;=±|k| i,j = l 



Re 



\ a J ' \ a a J \aJ \ a a 

\ a / ' \ a a / Va/ V a a/ 



1 \ r^ y^ 
1 + — I — — ^ Thb,|k| {r, y, t) 



t^ 






a" 



(4.18) 



|k| = 1 



The only substantial difference from the treatment of the 
middle-energy remainder lies in the prefactor generated 
by the param eteriz ation of the complex photon energy 
given in Eq. (4.17). The photon angular functions A^ 
and A^i^l {i,j = 1,2) for the high-energy part are defined 
in Eq. (5.8) of Ref. || and in Eq. (4.3) in JH for an 
arbitrary bound state. Special formulas for the ground 
state can be found in Eq. (5.9) of Ref. @. The func- 
tions Ak and A]^ are not identical to the photon angular 
functions for the infrared and middle-energy parts A^ 
{i,j = 1,2) which are used for the low-energy part of the 
calculation in Sec. HI. It might be worth mentioning that 
in p^-^, both the functions A^^ and A^^^ are denoted 
by the symbol A^ . It is clear from the context which of 
the functions is employed in each case. 



In the last line of Eq. ( 4.18| ), we implicitly define the 
terms Thb,|k[ as 

2 

TiiB,\t.\{r,y,t) = Y^ Y^ 

k;=±|k.| i:j = l 

-[/'(?)<4?'!-'")/4!) 

xA(^,^) 

\ a a/ 

With these definitions, the high-energy remainder can be 
rewritten as 

AShb = - / dt / dr / dySHB{r,y,t). (4.20) 

TT Jo Jo Jo 

There are no lower-order terms to subtract, and therefore 

A^HB = - ^-^ Fhb [nlj , Za) . 4.21 

For the high-energy remainder Fhb , the limits as Za — > 
read [see Eq. (4.15) in El] 



Fhb {nlj , Za) 



1 



21 + 1 



17 4 
In 2 Ui 

18 3 ' ■ 



^-21n2U 
2 ) K. 

5 \ n-2;- 1 
ln2 

6 In 



0{Za). 



(4.22) 



For the atomic states under investigation, this leads to 

-Fhb(1Si/2, Za) ^ -0.093 457 + 0{Za) , 

-Fhb(2Si/2, Za) = -0.023 364 + 0{Za) , 

Fhb(2Pi/2, Za) = 0.014 538 + 0{Za) , 

Fhb(2P3/2, Za) = -0.042 315 + 0{Za) . (4.23) 

The integration procedure for the high-energy part is 
adapted to the problem at hand. To this end, a crude 
estimate is found for the dependence of the function S'hb 
defined in Eq. (4.18) on its arguments. The consider- 



ations leading to this estimate are analogous to those 
outlined in Sec. HI A for the infrared part. The result is 



the approximate expression 

-y \ f^ ^ 
e " exp — 

for S'hb- This leads naturally to the definition 



(1 -r)y 



(4.24) 
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gHB = 1 



1 

at 



a-r), 



(4.25) 



so that the (approximate) dependence of S'hb on the ra- 
dial variable at large y is exp(— gHBj/)- Note that gnB 
may assume large values (^ 1) as t ^ 0; this is unlike 
the analogous quantity 



1 



(1 



We now turn to a brief discussion of the convergence 
acceleration techniques use d in t he evaluation of the func- 
tion S'hb defined in Eq. ( 4.18| ). The angular momen- 
tum decomposition of S'hb gives rise to a sum over the 
terms Thb,|k| [see the last line of Eq. (4.18)], where \k\ 
represents the modulus of the Dirac angular momentum 
quantum number of the virtual intermediate state. In 
shorthand notation, and suppressing the arguments, we 
have 



S'hb (X y^Ti 



HB,|k|- 



The radial Green function Gb = GB{ry/a,y/a, z) in co- 
ordinate space needs to be evaluated at the radial argu- 
ments r y / a and y / a (where < r < 1), and at the energy 
argument z = £■„ - w = i/2 {t~^ - t) [see Eq. ( [f.lSD ]. A 
crucial role is played by the ratio r of the two radial ar- 
guments. Indeed, for \k\ -^ oo, we have [see Eq. (4.7) 
inrai 



T, 



„2[k| 



HB,|« 



const. 



O 



(4.27) 



where "const." is independent of \n \ and depends only on 



r, y and t. The series in Eq. ( 4.26 ) is slowly convergent 
for r close to one, and the region near r = 1 is known to 
be problematic in numerical evaluations. Additionally, 
note that the region at r = 1 is more important at low 
Z than at high Z. This is because the function S'hb, for 
constant y, dep ends on r roughly as exp [—y (1 — r)/{a t)] 
[see Eq. (g3|)], where a = 2{Za)/n + 0[{Zaf]. For 
small Z, the Bohr radius l/{Za) of the hydrogenlike 
system is large compared to high-Z systems, which em- 
phasizes the region near r = 1. In this region, the se- 
ries in (4.26) is very slowly convergent. We have found 



that the convergence of this series near r = 1 can be ac- 
celerated very efhciently using the combined nonlinear- 
condensation transformation fUj^ app lied to the seri es 
Er=o *fe where tk = TBB,k+i [see Eqs. ( p^ ) and (pTJ) ]. 
We first transform this series into an alternating series 
by a condensation transformation due to Van Wijngaar- 
den rara, 



CXD OO 



k=o j=o 



in the infrared and the middle-energy part, where |c| = where 
|-\/l — z^l < 1 because z € (0,-E„). Having identified the 
leading exponential asymptotic behavior of the integrand 
S'hb, it is rather straightforward to evaluate the three- 
dimensional integral in Eq. ( |4jO| ) by Gauss-Laguerre and 
Gauss-Legendre q uad rature ||5|| [the scaling parameter a 
is defined in Eqs. (3.9) and ( p. 14 )]. The numerical results 
for the high-energy remainder function _Fhb are found in 
Table^. These results are consistent with the limits in 
Eq. (p|- 



T^^'h. 



0+i)-i- 



(4.28) 



(4.29) 



k=0 



We then accelerate the convergence of the alternating se- 
ries ^'^oi—^y Aj by applying the nonlinear delta trans- 
form (5„ (l,So), which is discussed extensively in [ p7[ . 
The explicit formula for this transformation is given by 
defining 



J=0 



(4.30) 



as the nth partial sum of the Van Wijngaarden trans- 
formed input series. The delta transform reads [see 
Eq. (8.4-4) of 1^]], 



(4,26) Sl^\l,So) 



EM) 

3=0 



n\ {l+j)n- 



jj (l-|-n)„_i Bj+1 



EM) 



i=o 



n\ (l+j)n- 



(4.31) 



jj (l + n)„_i Bj+i 



where 



B, = (-1)^A, 



(4.32) 



The convergence acceleration proceeds by calculating a 
sequence of transforms (S„ in increasing transformation 
order n. It is observed that the transforms con verge much 
faster than the partial su ms S„ defined in Eq. ( 4.30| ). The 
upper index zero in Eq. ( 4.31| ) indicates that the trans- 
formation is started with the first term Aq. 

The combined transformation (combination of the con- 
densation transformation and the Weniger transforma- 
tion) was found to be applicable to a wide range of slowly 
convergent monotone series (series whose terms have the 
same sign), and many examples for its application were 
given in Ref. llQ] . For the numerical treatment of radia- 
tive corrections in low-Z systems, the transformation has 
the advantage of removing the principal numerical diffi- 
culties associated with the slow convergence of angular 
momentum decompositions of the propagators near their 
singularity for equal radial arguments. 

In a typical case, sufficient pr ecisio n (10^^^) in the 
convergence of the sum in Eq. (4.26) is reached in a 
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transformation order n < 100 for the nonlinear trans- 
formation 6n (1, So), a region in which the nonhnear se- 
quence transformation S is numerically stable. Although 
the delta transformation exhibits considerable numeri- 
cal stability in higher transformation orders |lO,M, in- 
evitable round-off errors start to accumulate significantly 
in an excessively high transformation order of n « 500 
in a typical case |5[|, and this situation is avoided in the 
current evaluation because the transforms exhibit appar- 
ent convergence to the required accuracy before numer- 
ical round-off errors accumulate. Note that evaluation 
of the condensed series 



Aj in Eq. (4.29) entails sam- 



pling of terms ThbJkI for rather large \k\, while elimi- 
nating the necessity of evaluating all terms Tjjb.IkI up 
to the maximum index. The highest angular momentum 
\k\ encountered in the present calculation is in excess 
of 4 000 000. However, even in extreme cases less than 
3 000 evaluations of particular terms of the original se- 
ries are required. The computer time for the evaluation 
of the slowly convergent angular momentum expansion 
near the singularity is reduced by roughly three orders 
of magnitude by the use of the convergence acceleration 
methods. 

In certain parameter regions (e.g. for large energy of 
the virtual photon), a number of terms of the input series 
tk have to be skipped befo re th e convergence acceleration 
algorithm defined in Eqs. ( 4.28 )"( 4.32 ) can be appHed (in 
order to avoid transient behavior of the first few terms 
in the sum over k). In this case, the input data for the 
combined nonlinear-condensation transformation are the 
terms tk = Tub.Jc+i+k,, where Kg denotes the number of 
terms which are directly summed before the transforma- 
tion is applied. These issues and further details regarding 
the application of the convergence acceleration method to 
QED calculations can be found in Appendix H.2 of g]. 



C. Results for the High-Energy Part 

The limit of th e fun cti on Fa as Za -^ can be derived 
easily from Eqs. ( 4.14 ), ( 4.22 ) as a function of the bound 
state quantum numbers. For Fn the limit is 



Fii{nlj, Za) 
'11 



10 



-In 2 



16 
15 



- 2 In 2 



1 



2/ + 1 



3 1 „ 5 

In2 + - 

2 4 

17 4 , „ 
In 2 

18 3 



1 



k{21 + 1) 



3 , „ 9 

- In 2 - - 
2 4 



1 



OLO 



|-ln2 
6 



n-2l-l 
n{2l + l) 



-OiZa) 



(4.33) 



For the atomic states investigated here, this expression 
yields the numerical values 



Fh(1Si/2, Za) = -1.313 085 + 0{Za) , 
Fh(2Si/2, Za) = -1.446 418 + 0{Za) , 
i^H(2Pi/2, Za) = -1.066 667 + 0{Za) , 
i^H(2P3/2, Za) = -0.566 667 -I- 0{Za) . 
Numerical results for the high-energy part 
Fuinlj, Za) — F^Ainlj, Za) + FuBinlj, Za) 



(4.34) 



(4.35) 



are also summarized in Table [V. Note the apparent con- 



sistency of the numerical results in Table |IV| w ith their 
analytically obtained low-Z limits in Eq. ( [4.341) . 



TABLE V. Numerical results for the scaled self-energy function F and the self-energy remainder function Gse. 



z 


FilSr/2, Za) 


F{2Sr/2,Za) 


F{2P,^2,Za) 


F{2P3/2,Za) 


1 

2 
3 

4 
5 


10.316 793 659(1) 
8.528 325 061(1) 
7.504 503 432(1) 
6.792 824 089(1) 
6.251 627 086(1) 


10.546 825 185(5) 
8.758 870 25(1) 
7.735 777 20(1) 
7.025 002 41(1) 
6.484 860 42(1) 


-0.126 396 37(1) 
-0.125 816 16(1) 
-0.124 992 24(1) 
-0.123 968 79(1) 
-0.122 774 94(1) 


0.123 498 56(1) 
0.123 835 55(1) 
0.124 317 10(1) 
0.124 918 48(1) 
0.125 623 30(1) 


Z 


GsE(lSi/2,^a) 


GsE(2Si/2,^a) 


GsE(2Pi/2,^a) 


GsE(2P3/2,^a) 


1 
2 
3 

4 
5 


-30.290 24(2) 
-29.770 967(5) 
-29.299 169(2) 
-28.859 223(1) 
-28.443 372 3(8)'' 


-31.185 15(9) 
-30.644 66(5) 
-30.151 93(2) 
-29.691 27(1) 
-29.255 033(8) 


-0.973 5(2) 
-0.949 40(5) 
-0.926 37(2) 
-0.904 12(1) 
-0.882 478(8) 


-0.486 5(2) 
-0.470 94(5) 
-0.456 65(2) 
-0.443 13(1) 
-0.430 244(8) 



''The result for this entry given in W contains a typographical error. 
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V. COMPARISON TO ANALYTIC 
CALCULATIONS 

The numerical results for th e sc aled self-energy func- 
tion F{nlj,Za) defined in Eq. (2.1) are given in Table 0, 
together with the results for the nonperturbative self- 
energy remainder functi on GsEJnli, Za), which is im- 
plicitly defined in Eq. (2.2). Results are provided for 
K- and L-shell states. The results here at Z = 5 are 
consistent with and much more precise than the best 
previous calculation |Q. The numerical results for the 
self-energy remainder Gse are obtained by subtracting 
the analytic lower-order terms listed in Eq. (2^) from 
the complete numerical result for the scaled self-energy 
function F(nlj, Za). No additional fitting is performed. 

Analytic and numerical results at low Z can be com- 
pared by considering the self-energy remainder function 
Gse- Note that an inconsistency in any of the ana- 
lytically obtained lower-order terms would be likely to 
manifest itself in a grossly inconsistent dependence of 
GsEinlj, Za) on its argument Za; this is not observed. 
For S states, the following analytic model for Gse is 
commonly assumed, which is motivated in part by a 
renormalization-group analysis |58|] and is constructed in 
analogy wi th the pat tern of the analytic coefficients Aij 
in Eq. (U) and (pj) 

GsEinSi/2, Za) = ^60(^81/2) 

+ iZa) [A7i(nSi/2) HZa)-^ + AninSi/2)] 

+ (Za)2 [As3{nSi/2) ln\Za)'^ 

+A82(nSi/2) In^iZa)-^ 

+Asi{nSif2) HZa)-^ + AsoinSi/2)] . (5.1) 

The (probably nonvanishing) A^s coefficient, which intro- 
duces a triple logarithmic singularity at Za = 0, hinders 
an accurate comparison of numerical and analytic data 
for Gse- A somewhat less singular behavior is expected 
of the difference 

AGsE(^a) = Gse(2Si/2, Za) - Gse(1Si/2, Za) , (5.2) 

because the leading logarithmic coefficients in any given 
order of Za are generally assumed to be equal for all S 
states, which would mean in particular 



A7l(lS 



I/2J 



A7i(2Si/2) and 



^83(181/2) = ^83(281/2) 



(5.3) 



Now we define AAki as the difference of the values of the 
analytic coefficients for the two lowest S states: 



AAki^Akii2S 



1/2J 



AkiilS 



1/2J 



(5.4) 



The function AGse defined in Eq. ( |5.2D can be as- 
sumed to have the following semi-analytic expansion 
about Za = 0: 



AGsEiZa) = AAeo + (Za) AA70 
+ {Za)^ [AAs2 \n^{Za)-^ 

+AAfii ln{Za)-^ + AAsa + o{Za)] 



(5.5) 



In order to detect possible inconsistencies in the numer- 
ical and analytic data for Gse, we difference the data 
for AGsE, i-G-, we consider the following finite difference 
approximation to the derivative of the function AGse^ 



g{Z) = AGseHZ + l)a)- AGse(Z«) ■ 



(5.6) 



We denote the analytic and numerical limits of 
AGsE{Za) as Za ^ as AA\^^' and AAg™ , respec- 
tively, and leave open the possibility of an inconsistency 
between numerical and analytic data by keeping AAg™ 

and A^gp" as distinct variables. In order to illustrate 
how a discrepancy could be detected by investigating the 
function g{Z), we consider special cases of the function 
AGsE{Za) and g{Z). We have for Z = 0, which is de- 
termined exclusively by analytic results, 



AGse(O) = AA 



(an) 
60 



(5.7) 



whereas for Z — 1, which is determined by numerical 
data. 



AA 



(„u) 



60 



a [AA' 



70 



Ola) 



(5.8) 



AGsE(a) 
and for Z — 2, 

AGsE(2a) = a4"o") + a [2 AAro + o(a)] , (5.9) 
etc. Hence for Z = 0, we have 
g(0) = AGsE(a) - AGse(O) 



— A /)('!") 



AA 



60 



AA. 



(an) 
60 



\AA' 



70 



o{Za)] 



(5.10) 



For Z ~ 1, the value of g is determined solely by numer- 
ical data, 

ff(l) = AGsE(2a)-AGsE(a) 

= a [AAjo + oiZa)] , (5.11) 

and for Z — 2, we have 

ff(2) = AGsE(3a)~AGsE(2a) 

= a [AAjo + o{Za)] . (5.12) 

Analogous equations hold for Z > 2. The analytic and 
the numerical data from Table M lead to the five values 
5(0), g{l), 5(2), g(3), and g(4). A plot of the function 
g{Z) serves two purposes: First, the values 5(1), . . . , g(4) 
should exhibit apparent convergence to some limiting 
value a AA-^q as Z ^ 0, and this can be verified by in- 
spection of the plot. Secondly, a discrepancy between the 
analytic and numerical approaches would result in a non- 



(nu) 
60 



AAgg" which would appear 



vanishing value for AA, 

as an inconsistency between the trend in the values of 

5(1), . . . , and 5(4) and the value of g{0) [see Eq. (5.10)]. 
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FIG. 4. Plot of the function g{Z) defined in Eq. (|5J) in 
the region of low nuclear charge. For the evaluation of the 
data point at Z = 0, a value of ^60(181/2) = -30.924 15(1) 
is employed @,|l5[||. 
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FIG. 5. Comparison of numerical data and analyti- 

cally evaluated higher-order binding corrections for t he 2P 1/2 
state. We plot the function (71/2 (■^) defined in Eq. (5.15) in 
the region of low Z. The numerical data obtained in the cur- 
rent investigation appear to be consistent with the analytic 
result of y4Bo(2Pi/2) = -0.998 91(1) obtained in |l6). 



Among the separate evaluations of Aqq for the ground 
state, the result in |]l5| has the smallest quoted uncer- 
tainty. In Fig. ^ we display a plot of g{Z) for low nu- 
clear charge Z. A value of ^60(181/2) = Aq^'{1'&i/2) 



-30.92415(1) i,[ll|9| is used in Fig. |. The results in- 
dicate very good agreement between the numerical and 
analytic approaches to the Lamb shift in the low-Z re- 
gion up to the level of a few Hz in frequency units for the 
low-lying atomic states (where n is the principal quan- 
tum number). The error bars represent the numerical 
uncertainty of the values in Table M, which correspond 
to an uncertainty on the level of 1.0 x Z'^ Hz in frequency 
units. 



Analytic work on the correction Ago has extended over 
three decades |15|, p4|-p7| . The complication arises that 
although the calculations are in general analytic, some re- 
maining one-dimensional integrations could not be evalu- 
ated analytically because of the nature of the integrands 
[see e.g. Eq. (6.96) in |15)]. Therefore a step-by-step com- 
parison of the analytic calculations is difficult. An addi- 
tional difficulty is the isolation of those analytic terms 
which contribute in a given order in Za, i.e., the isola- 
tion of only those terms which contribute to Ago- The 
apparent consistency of the numerical and analytic data 
in Fig. H represents an independent consistency check on 
the rather involved analytic calculations. 



Our numerical results are not inconsistent with the an- 
alytic result Ipl for a higher-order logarithm, 



A 



71 



139 
"64" 



In 2 



4.65, 



(5.13) 



although they do not necessarily confirm it. As in [H, 
we obtain as an estimate ^471 = 5.5(1.0) (from the fit to 
the numerical data for both S states). Logarithmic terms 
corresponding to the (probably) nonvanishing Ag,^ coef- 
ficient should be taken into account for a consistent fit 
of the corrections to Gse- These highly singular terms 
are difficult to handle with a numerical fitting procedure. 
The terms Aga, A^2 and Agi furnish three more free pa- 
rameters for the numerical fit, where only five data points 
are available (in addition to the quantities AgOi ^71 and 
Ayo, which may also be regarded as free parameters for 
the fitting procedure). The determination of Aqq by a fit 
from the numerical data is much more stable than the de- 
termination of the logarithmic correction A^i . We briefly 
note that our all-order evaluation essentially eliminates 
the uncertainty due to the unknown higher-order ana- 
lytic terms. Also, it is interesting to note that the same 
numerical methods are employed for both the S and P 
states in our all-order (in Za) calculation, whereas the 
analytic treatment of S and P states differs [|l5|, Q . 

The comparison of numerical and analytic results is 
much less problematic for P states, because the function 
Gse is less singular [see Eqs. (2^) and (|^)]. For the 2P 
states, we observe that the function GsE,{'2Pj,Za) has 
the same semi-analytic expansion about Za = as the 
function /S.Gs^{Za) defined for S states in Eq. (^.2|). We 
have 
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GsE{2Fj,Za) = Aeo{'2Fj) + iZa)Ajo{2Pj) 

+ {Zaf [v482(2Pj) IniZay 

+A8i(2P,) ln(Za)-2 + A8o(2Pj) + o(Za)] . (5.14) 
Hence, we plot the function 

g,{Z) = Gse(2P,, {Z + l)a)- Gse(2Pj, Za) (5.15) 

for j = 1/2 and j = 3/2 in the region of low Z, with the 
notion that an inconsistent analytic result for yl6o(2Pj) 
would lead to irregularity at Z = 0, in analogy with the 
S states. The numerical data shown in Figs. ^, and g 
appear to be consistent with the analytic results of 



^6o(2Pi/2) = -0.998 91(1) and 



Aao(2P 



zji) 



-0.503 37(1) 



(5.16) 



obtained in (l6|. In this context it may be interesting 
to note that analytic results obtained in [|6|^ for the 
higher-order binding corrections to 2P, 3P, and 4P states 
have recently been confirmed indirectly [pffl. Finally, al- 
though it may be possible to obtain more accurate esti- 
mates of some higher-order analytic corrections, notably 
the AiQ coefficient for P states and AAyg for the two 
lowest-lying S states, we have not made such an analysis 
in the current work; we have restricted the discussion to 
a check of the consistency with the available results for 

^60- 
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FIG. 6. For the 2 P3/2 state, we plot the function g^/2{Z) 
defined in Eq. (5.15) in the region of low Z. The numerical 



data obtained in the current investigation appear to be con- 
sistent with the analytic result of ^6o(2P3/2) = —0.503 37(1) 
from ^. 



two-body corrections at low Z and the comparison of an- 
alytic and numerical results [p8| , |6l| -i72| . Traditionally, the 
self-energy correction for hydrogenlike systems has posed 
a computational challenge. Here we have described a 
nonperturbative evaluation of the one-photon self-energy 
correction in hydrogenlike ions with low nuclear charge 
numbers Z — 1 to 5. The general outline of our approach 
is discussed in Sec. ||. In Sec. HI, the numerical evalua- 
tion of the low-energy part (generated by virtual photons 
of low energy) is described. In Sec. IV, we discuss the nu- 
merical evaluation of the high-energy part, which is gen- 
erated by high-energy virtual photons and contains the 
formally infinite contributions, which are removed by the 
renormalization. Sec. IV also contains a brief discussion 



of the convergence acceleration methods as employed in 
the current evaluation. We discuss in Sec. M the com- 
parison of analytic and numerical data for K- and L-shell 
states in the region of low Z. The main results of this 
paper are contained in Table |y|: numerical data, non- 
perturbative in Za, for the scaled self-energy function F 
and the self-energy remainder function Gse for K- and 
L-shell states at low nuclear charge. The numerical ac- 
curacy of our data is 1 Hz or better in frequency units 
for IS, 2S and both 2P states in atomic hydrogen. 

The comparison of analytic and numerical results to 
the level of accuracy of the numerical data, which is dis- 
cussed in Sec. ^ indicates that there is very good agree- 
ment for the K- and L-shell states. The analytic and 
numerical data are shown in Figs. 0, 0, and O. Our all- 
order evaluation eliminates any uncertainty due to the 
unknown higher-order analytic terms; the current numer- 
ical uncertainty in the self energy is at the level of 1 Hz 
for atomic hydrogen. 
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VI. CONCLUSION 

There has recently been a rather broad interest in the 
numerical calculation of relativistic, QED self energy, and 
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